Diversity Tackled with R
Look at your fingers; controlled by the mind can do great things. However, imagine if each one has a little brain of its own, with different ideas, desires, and fears ¡How wonderful things will be made out of an artist with such hands! -Ode to multidisciplinarity
First plunge into diversity
Species diversity, in its simplest definition, is the number of species in a particular area and their relative abundance (evenness).
Once we know the taxonomic composition of our metagenomes, we can do diversity analyses. Here we will discuss the two most used diversity metrics, α diversity (within one metagenome) and β (across metagenomes).
- α Diversity: Can be represented only as richness (, i.e., the number of different species in an environment), or it can be measured considering the abundance of the species in the environment as well (i.e., the number of individuals of each species inside the environment). To measure α-diversity, we use indexes such as Shannon’s, Simpson’s, Chao1, etc.
- β diversity is the difference (measured as distance) between two or more environments. It can be measured with metrics like Bray-Curtis dissimilarity, Jaccard distance, or UniFrac distance, to name a few. Each one of this measures are focused on a characteristic of the community (e.g., Unifrac distance measures the phylogenetic relationship between the species of the community).
In the next example, we will look at the α and the β components of the diversity of a dataset of fishes in three lakes. The most simple way to calculate the β-diversity is to calculate the distinct species between two lakes (sites). Let us take as an example the diversity between Lake A and Lake B. The number of species in Lake A is 3. To this quantity, we will subtract the number of these species that are shared with the Lake B: 2. So the number of unique species in Lake A compared to Lake B is (3-2) = 1. To this number, we will sum the result of the same operations but now take Lake B as our reference site. In the end, the β diversity between Lake A and Lake B is (3-2) + (3-2) = 2. This process can be repeated, taking each pair of lakes as the focused sites.
If you want to read more about diversity, we recommend to you this paper on the concept of diversity.
α diversity
β diversity
Diversity β measures how different two or more communities are, either in their composition (richness) or in the abundance of the organisms that compose it (abundance).
- Bray-Curtis dissimilarity: The difference in richness and abundance across environments (samples). Weight on abundance. Measures the differences from 0 (equal communities) to 1 (different communities)
- Jaccard distance: Based on the presence/absence of species (diversity). It goes from 0 (same species in the community) to 1 (no species in common)
- UniFrac: Measures the phylogenetic distance; how alike the trees in each community are. There are two types, without weights (diversity) and with weights (diversity and abundance)
There are different ways to plot and show the results of such analysis. Among others, PCA, PCoA, or NMDS analysis are widely used.
In the next picture, there are two lakes with different fish species:
Which of the options below is true for the alpha diversity in lakes A, B, and beta diversity between lakes A and B, respectively?
- 4, 3, 1
- 4, 3, 5
- 9, 7, 16
Please, paste your result on the collaborative document provided by instructors.
Plot alpha diversity
We want to know the bacterial diversity, so we will prune all non-bacterial organisms in our merged_metagenomes Phyloseq object. To do this, we will make a subset of all bacterial groups and save them.
> merged_metagenomes <- subset_taxa(merged_metagenomes, Kingdom == "Bacteria")Now let us look at some statistics of our metagenomes. By the output of the sample_sums() command:
> sample_sums(merged_metagenomes)
0-kraken_report 1-kraken_report 2-kraken_report
32389 53156 55956we can see how many reads there are in the library.
Library B_Sample_97 is the smallest with 32389 reads, while library B_Sample_99 is the largest with 55958 reads.
Also we can obtain the Max, Min, and Mean output on summary(), which can give us a sense of the evenness. For example, the OTU that occurs the most in the sample B_Sample_97 occurs 5979 times, while on average in sample B_Sample_98, an OTU occurs in 48.02 reads.
summary(merged_metagenomes@otu_table@.Data)
0-kraken_report 1-kraken_report 2-kraken_report
Min. : 0.00 Min. : 0.00 Min. : 0.00
1st Qu.: 0.00 1st Qu.: 1.00 1st Qu.: 1.00
Median : 1.00 Median : 2.00 Median : 2.00
Mean : 29.26 Mean : 48.02 Mean : 50.55
3rd Qu.: 3.00 3rd Qu.: 7.00 3rd Qu.: 6.50
Max. :5979.00 Max. :4459.00 Max. :7048.00 To have a more visual representation of the diversity inside the samples (i.e., α diversity), we can now look at a graph created using Phyloseq:
> plot_richness(physeq = merged_metagenomes,
measures = c("Observed","Chao1","Shannon")) 
Each of these metrics can give an insight into the distribution of the OTUs inside our samples. For example, the Chao1 diversity index gives more weight to singletons and doubletons observed in our samples, while Shannon is an entropy index remarking the impossibility of taking two reads out of the metagenome “bag” and that these two will belong to the same OTU.
The code and the plot using the three options will look as follows:
The “title” option adds a title to the figure.
> plot_richness(physeq = merged_metagenomes, title = "Alpha diversity indexes for three samples from Healthy Upper Epidermidis", measures = c("Observed","Chao1","Shannon"))
The “nrow” option arranges the graphics horizontally.
> plot_richness(physeq = merged_metagenomes, measures = c("Observed","Chao1","Shannon"), nrow=3) 
The “sortby” option orders the samples from least to greatest diversity depending on the parameter. In this case, it is ordered by “Shannon” and tells us that the B_Sample_97 has the lowest diversity and the B_Sample_98 the highest.
> plot_richness(physeq = merged_metagenomes, measures = c("Observed","Chao1","Shannon"), sortby = "Shannon") 
Considering those mentioned above, together with the three graphs, we can say that B_Sample_98 and B_Sample_99 present a higher diversity compared to sample B_Sample_97. Moreover, the diversity of the sample B_Sample_98 and B_Sample_99 is mainly given by singletons or doubletons. While the diversity of sample B_Sample_98 provided by species occurs in greater abundance. Although the values of H (Shannon) above three are considered to have a lot of diversity.
A caution when comparing samples is that differences in some alpha indexes may be the consequence of the difference in the total number of reads of the samples. A sample with more reads is more likely to have more different OTUs, so some normalization is needed. Here we will work with relative abundances, but other approaches could help reduce this bias.
Absolute and relative abundances
From the read counts that we just saw, it is evident that there is a great difference in the number of total sequenced reads in each sample. Before we further process our data we should look to see if we have any non-identified reads. Marked as blank (i.e.,““) on the different taxonomic levels:
> summary(merged_metagenomes@tax_table@.Data== "")
Kingdom Phylum Class Order
Mode :logical Mode :logical Mode :logical Mode :logical
FALSE:1107 FALSE:1107 FALSE:1107 FALSE:1106
TRUE :1
Family Genus Species
Mode :logical Mode :logical Mode:logical
FALSE:1024 FALSE:860 TRUE:1107
TRUE :83 TRUE :247
With the command above, we can see blanks on different taxonomic levels. For example, we have 247 blanks at the genus level. Although we could expect to see some blanks at the species or even at the genus level; we will get rid of the ones at the genus level to proceed with the analysis:
> merged_metagenomes <- subset_taxa(merged_metagenomes, Genus != "") #Only genus that are no blank
> summary(merged_metagenomes@tax_table@.Data== "")
Kingdom Phylum Class Order
Mode :logical Mode :logical Mode :logical Mode :logical
FALSE:860 FALSE:860 FALSE:860 FALSE:859
TRUE :1
Family Genus Species
Mode :logical Mode :logical Mode:logical
FALSE:856 FALSE:860 TRUE:860
TRUE :4 Next, since our metagenomes have different sizes, it is imperative to convert the number of assigned reads (i.e., absolute abundance) into percentages (i.e., relative abundances) to compare them.
Right now, our OTU table looks like this:
> head(merged_metagenomes@otu_table@.Data)
0-kraken_report 1-kraken_report 2-kraken_report
46157 5979 4459 7048
2444 236 308 517
26042 94 98 207
26048 10 40 34
26054 7 19 17
26040 7 25 29To make this transformation to percentages, we will take advantage of a function of Phyloseq:
> percentages <- transform_sample_counts(merged_metagenomes, function(x) x*100 / sum(x) )
> head(percentages@otu_table@.Data)
head(percentages@otu_table@.Data)
0-kraken_report 1-kraken_report 2-kraken_report
46157 23.78092435 12.09384323 17.36559405
2444 0.93866836 0.83536751 1.27383827
26042 0.37387638 0.26579875 0.51002809
26048 0.03977408 0.10848929 0.08377273
26054 0.02784186 0.05153241 0.04188636
26040 0.02784186 0.06780580 0.07145321Now, we are ready to compare the abundaces given by percantages of the samples with beta diversity indexes.
Beta diversity
As we mentioned before, the beta diversity is a measure of how alike or different our samples are (overlap between discretely defined sets of species or operational taxonomic units). To measure this, we need to calculate an index that suits the objectives of our research. By the next code, we can display all the possible distance metrics that Phyloseq can use:
> distanceMethodList
$UniFrac
[1] "unifrac" "wunifrac"
$DPCoA
[1] "dpcoa"
$JSD
[1] "jsd"
$vegdist
[1] "manhattan" "euclidean" "canberra" "bray" "kulczynski"
[6] "jaccard" "gower" "altGower" "morisita" "horn"
[11] "mountford" "raup" "binomial" "chao" "cao"
$betadiver
[1] "w" "-1" "c" "wb" "r" "I" "e" "t" "me" "j" "sor" "m"
[13] "-2" "co" "cc" "g" "-3" "l" "19" "hk" "rlb" "sim" "gl" "z"
$dist
[1] "maximum" "binary" "minkowski"
$designdist
[1] "ANY"
Describing all these possible distance metrics is beyond the scope of this lesson, but here we show which are the ones that need a phylogenetic relationship between the species-OTUs present in our samples:
- Unifrac
- Weight-Unifrac
- DPCoA
We do not have a phylogenetic tree or phylogenetic relationships. So we can not use any of those three.
We will use Bray-curtis since it is one of the most robust and widely used distance metrics to calculate beta diversity.
Let’s keep this up! We already have all we need to begin the beta diversity analysis. We will use the Phyloseq command ordinate to generate a new object where the distances between our samples will be allocated after calculating them. For this command, we need to specify which method we will use to generate a matrix. In this example, we will use Non-Metric Multidimensional Scaling or NMDS. NMDS attempts to represent the pairwise dissimilarity between objects in a low-dimensional space, in this case, a two-dimensional plot.
meta_ord <- ordinate(physeq = percentages, method = "NMDS", distance = "bray")If you get some warning messages after running this script, fear not. It is because we only have three samples. Few samples make the algorithm warn about the lack of difficulty in generating the distance matrix.
By now, we just need the command plot_ordination() to see the results from our beta diversity analysis:
> plot_ordination(physeq = percentages, ordination = meta_ord)
In this NMDS plot, each point represents the combined abundance of all its OTUs. As depicted, each sample occupies space in the plot without forming any clusters. This output is because each sample is different enough to be considered its own point in the NMDS space.